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1.  Introduction 


lat  ua  consider  the  problen  of  finding  a  fwcilm  *U»y)  which  satiaflaa 
th*  parabolic  partial  differential  e»juatioo 

(1.1) 

in  Ri  0<x  <k,  0  <  t  <T  and  which  satiafiaa  tba  Initial  cooditicna 

(1.2)  uU,0)  ■'«(*) 

and  the  boundary  ccnditicna 

a^tJu^O.t)  ♦  ^(OotO.t) 

OgttJu^CAjt)  ♦  P2(t)uU,t) 

Tha  function  u(x,y)  la  r^uirad  to  bo  ccntlnaoua  and  bowdad  in  R  ♦  S,  whara  8 
i,  the  boundary  of  R,  except  that  continuity  of  »Ujr)  i*  not  ra^uirad  at  a 
point  of  S  where  a  Junp  diacontinuity  occura  in  any  of  tha  other  functions  i»- 
volred  in  (1.2)  oral  (1.5).  *e  assuns  that  aU  functions  ara  continuoua  arc 

bounded  in  R. 

in  analytic  solution  for  tha  abort  problem  ia  arailabla  only  undar  wary 
special  conditions,  ta  shall  not  sssk  aooh  solutions  hart.  Koraoaar,  wa  shall 
not  even  consider  the  rary  laportant  nutations  of  aaiatanoa  and  unWuenaas. 
Bather,  we  stall  be  concerned  with  the  study  of  numerical  procedures  for  aolw- 
tnc  certain  difference  ^uation.  which  one  obtain,  ona  uaa.  finite  differ¬ 
ence  mthoda.  Of  course  it  is  recognised  that  ary  eolations  ao  cfcUlnad  would 
ba  of  no  Talus  in  tin  erant  that  no  solution  of  tha  criminal  prcblan  axiatad 
or  in  tha  treat  Uni  nars  than  one  solution  axiatad.  In  certain  cans*  axiatano. 


0<  x 

-  $(t)  0  *  t 

-  f(t)  0*t. 


uid  uniqueness  ara  available  (See  for  iwtauM  III  and  [ft]  ). 

The  numerical  methods  which  —  shell  on— ider  lndnda  tbc  well-knoua 

forward  Difference  method,  the  Cnrik-ki  nolson  — thod  £1*1  and  the  Dufc’ort-JFrankd 
method  [51*  With  the  Crank-hicoleon  —thod  one  is  usually  required  to  moire 

at  each  step  a  system  of  simultaneous  equations.  This  can  be  done  effectively 

in  mar  cases  using  the  Successive  Over-relaxation  —thod  of  iteration  [22]. 

Alternatively  if  the  ays  tee  la  linear,  or  ia  replaced  by  a  related  linear 

system,  one  can  use  a  non-iterative  method,  apparently  first  used  for  this 

type  of  problem  by  ft-uce,  Pe  iceman,  Rachfard,  and  Rice  [2]. 

In  section  2  we  first  derive  finite  difference  representations  of  the 
boundary  conditions  including  two  well-known  methods  and  one  which  we  have 
not  seen  in  the  literature.  Ve  than  derive  several  difference  equations  end 
obtain  local  error  estimates  using  methods  of  numerical  quadrature.  Kext,  we 
consider  the  process  known  as  "extrapolation  to  sero  grid  sIm",  fit's t  used 
by  L.F.  Richardson  [20j.  Using  this  process  one  can,  under  faro rabbi  condi¬ 
tions,  appreciably  Lap rove  the  accuracy  of  the  finite  difference  nothods. 

We  indicate  how  results  obtained  using  several  different  grid  si  mi  can  often 
be  used  to  obtain  a  more  accurate  extrapolation  formula  than  one  would  obtain 
using  Richardson's  forwla  wherein  it  is  essi— d  that  the  error  tends  to  sero 
like  the  square  of  the  grid  else. 

In  order  to  cotqure  the  varlo—  finite  difference  amtbods  we  prepared 
a  comprehensive  ratine  for  Ordvac  which  is  capable  of  solving  a  large  class 
of  problems  by  each  —thod.  Details  of  this  program  are  given  in  section  l» 
end  Appendix.  The  apedflx  problems  which  —  considered  are  discussed  in 
section  3« 
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Tha  progr**  r>»sy  be  considered  u  a  research  tool  to  enable  on*  to  study 
the  effectiveness  of  the  various  Methods  under  different  conditions.  When  need 
together  with  a  sore  ooaplete  theoretical  analysis  than  we  has*  been  able  to 
carry  out  here,  it  is  to  be  expected  that  new  and  significant  results  co 
finite  difference  Methods  should  be  obtained. 

In  section  5  the  results  obtained  on  the  Ordvae  are  analysed  and  tentative 
conclusions  axe  reached  as  to  the  order  of  convergence  and  the  overall  effective¬ 
ness  of  each  of  the  net hods.  For  the  linear  case  our  results  indicate  strongly 
that  the  use  of  the  Crank-Nloolson  method  with  the  non-iterative  procedure, 
is  superior  to  any  of  the  other  aetbods.  In  the  non-linear  case  however,  one 
loses  accuracy  with  this  Method  if  one  Modifies  the  difference  equation  ir  an 
attempt  to  avoid  all  iterations.  Further  work  appears  to  be  indicated  in 
order  to  preserve  the  accuracy  of  the  difference  equation  without  requiring 
too  many  iterations.  Vie  suggest  one  possible  schme  for  doing  this,  but  we 
have  not  yet  tried  it  out  on  the  Machine. 

For  the  Forward  Difference  Method  and  for  the  C rank -Nlco Ison  Method,  used 
with  the  iterative  process,  we  found  that  the  accuracy  could  usually  be  In- 
pro  >wd  considerably  by  the  use  of  tha  proper  extrapolation  to  sero  grid  else. 
Such  extrapolation  was  more  difficult  to  carry  out  for  the  Crank-Klcolscai 
method, used  without  iterations. 

As  expected  the  use  of  the  proposed  Method  fbr  representing  the  boundary 
conditions  appeared  to  be  More  accurate  than  either  of  the  other  two  considered. 

Although  this  study  is  far  f  nan  complete,  it  was  felt  that  the  work 
which  has  been  done  to  date  should  be  reported  now,  especially  since  no  fur¬ 
ther  work  can  be  done  daring  the  next  acadwAc  year.  Tha  work  may  be  regarded 
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as  a  natural  extension  of  research  which  was  begun  tgr  the  Interior  Ballistics 
Laboratory,  Aberdeen  Pro  ring  Ground.  One  «f  the  objects  m a  to  determine  the 
order  of  convergence  of  the  forward  Difference  Method.  Extensive  calculations 
were  carried  out  on  the  Bell  Rolay  Corputer,  as  daecrlbed  in  [6],  and  on  the 
Eniac.  The  latter  calculations  have  not  previously  been  reported.  We  shell 
study  then  in  section  5  along  with  the  Ordvac  results. 

In  addition  to  acknowledging  the  sponsorship  of  the  Office  of  Ordnance 
Research,  U.S.  Amy,  we  should  also  like  to  express  our  appreciation  for 
tbs  encourage  went  and  assistance  fumivhed  by  the  Interior  Ballistic  Labora¬ 
tory  at  the  Aberdeen  Proving  Ground.  Discussions  with  Dr.  V.C.  Taylor  of  the 
Interior  Ballistic  Laboratory  end  with  Dr.  B.L.  Hides,  forwerly  of  that  lab¬ 
oratory,  were  of  great  value.  We  are  also  pleased  to  acknowledge  the  assis¬ 
tance  of  Dr.  R.M.  Coni  an,  Hiss  J.M.  Wood,  Ur.  B.T.C.  Koo,  and  Ur.  3.  Soscia 
who  worked  at  various  tises  on  the  prograrr.ij.^  :f  the  problea  for  Ordvac. 


2.  Difference  aquations 


Ut  h.  and  k  be  positive  ousters  such  that  */h  K  an  integer,  which  ee 
denote  by  1,  and  let  denote  the  act  of  point*  (x,t)  such  that  x/h  and 
t/k  are  Integer*.  let  R^*  and  denote  respectively  the  eet  of  points 

of  1^  k  belonging  to  R  and  to  3, (See  figure  1).  !e  replace  the  problem  of 
the  previous  section  by. that  of  finding  a  function  u(x,t)  defined  on  Rj^jj  *  ®h,k 
satisfying  a  certain  difference  equation  on  R^  ^  and  satisfying  on 
certain  otter  condition*  derived  froa  (1.2}  and  (1.3) • 


In  thi*  section  ee  derive  finite  difference^  which  '-e  obtained  froa  the 
differential  equation  and  froa  the  boundary  conditions,  and  »eek  to  ratlnato 
the  accuracy.  Our  aethoda  are  by  no  a* ana  rigorous-  Indeed  they  are  aerely 
foraal  aethoda  utlliaing  Taylor  aeries  expan* loo*  ehenever  convenient,  never¬ 
theless,  It  is  felt  that  these  aethoda  will  give  soar  basis  for  Baking  com¬ 
parison*  of  tte  ac  curse  lei  of  th*  various  finite  differeno*  aethod^at  least 


for  sqm  cases.  «e  will  later  tsat  ths  comparisons  thus  obtained  against 
numerical  results  obtained  for  typical  problsM  on  the  Ordme. 

For  discussion  of  finite  difference  analogues  of  problsM  involving 
infinite  region  —<«<«*>•  the  reader  is  referred  to  papers  by  tiourant, Fried¬ 
richs  -nd  Lewy  p]  and  by  fits  John  p<| . 

In  addition  to  deriving  the  difference  stations  we  will  alao  describe 
here  the  associated  naeerical  procedures.  In  order  that  we  Mr  be  able  to 
glre  a  c caplets  description  of  these  procedures  it  is  convenient  to  firwt 
consider  the  treatment  of  the  bounder}-  conditions. 
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2.1  Bbundary  condition*.  Although  our  discussion  could  bo  extended 


to  cover  no re 

general  situations,  we  shall  conaidsr 

only  tee  following  tee 

types  of  boundary  conditions 

_ 

r(.) 

u!o,t)  •  #t) 

,  (t  JO) 

(2.1.1) 

((b) 

,  (tso). 

u(A,t)  •  ftt) 

and 

(2.1.2) 

fu> 

u  (0,t)  *  -H(u 

-  u(0,t) )  , 

(t»o) 

j 

x  g 

( (b) 

ux(A,t)  -  0 

$ 

(t»0)  . 

Here  H  and  u 

ere  constants. 

be  shall  at  tinss  refer  to  the  boundary  conditions  represented  by  (2.1.1)  and 

I 

1 

(7.1.2)  as  B  V  and  N  D  conditions,  respectively. 

■*  first  note  that  the  Initial  conditions  (1.2)  are  easily  represented  by 

(2.1.3)  a,  •  g4  1  *  i  *M-1  . 

Here  nnd  subsequently  we  let  uif  jMih, »*  K±mt (i*),ete.,  teere  i  and  J  ere 
Integers.  Moreover,  the  B  7  conditions  (2.1.1)  can  be  replaced  by 

s 

(2.1.ltj  j  -  Pi  ,  (0  *  J) 

j  \\im*s  •  t0‘J)  • 

I  However,  with  the  HE  condition  (2.1.2a)  we  shall  consider  several  alter- 

j  native  procedures.  A  nunber  of  these  ars  c  on  aide  red  by  Prloe  and  Slack  £93. 

v«  shill  discuss  t*o  of  t!u»c«  to^etlier  with  one  which  nay  very  vail  be  known 
also,  although  we  have  not  seen  it  in  the  literature, 
i  In  order  to  simplify  our  notation  we  (hall  consider  a  configuration  of 

four  points  mar  the  line  x  *  0,  as  shown  in  ligure  2«  be  denote  by  e^, 

I 

•  x4,  f.  the  values  of  the  respective  variables  at  the  point  cuteered  1,  (1  •  0, 

i  *  * 

|  1,  4,  S).  be  seek  relation*  between  u^  and  one  or  sore  at  the  ether  e^.  It 
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differential  station.  Instead  we  shall  atteept  to  infer  this  froa  the  results 
of  Ordvac  c  Deputations. 

If  une  nukes  use  of  the  differential  equation  at  the  point  (x4,t4) ,  one 
can  obtain  wore  accurate  f  oracle  a.  Thus,  using  a  Taylor  caries  expansion 
about  (x4,t4)  we  bare 

*b  *  \  4  k"t 4  *k%t  4*“ 

Ug  •  u4  ♦  ta  ♦  Jh^  ♦  eh3  xSDat  ♦  •  •  • 

Here  the  derivatives  are  evaluated  at  (x4,t4).  Using  (1.1)  and  the  relation 
u  -  -  H(u  -  u i  )  we  obtain 

*  c  * 

(2.1.7)  t/i  w^d-Zr)  ♦  2nig  ♦  2rhH(u^-  u4>  ♦  kf4 

where  tfcrt  leading  terse  of  U.e  difference  between  the  rlglA  and  the  left  set¬ 
ter;  of  (2.1.7) ,  that  is,  the  ri^ht  laenLer  sdnus  tie  left  seaber,  are 

(!.!.»)  -T*«.  *  £«tf 

here  and  subsequently  *•  let 

(2.1.9)  r  -  \  . 

b* 

Price  and  SI.  ck  [lsj  show  tlat  thie  fonsila  introduces  instability  for  r  >  J, 
in  contrast  to  (2.1.5)  which  is  stable  for  al)  r.  Titus  tssusing  k  •  rh  ,  where 
r  1c  a  constant  not  greater  tUn  },  we  obtain  froa  (2.1.6) 

(2.1.1D)  -^'xxx  *  * 

Thus  the  leading  tars  of  tie  local  error  is  0(h3)  as  compared  with  0(h2)  with 
Option  1.  thforturwtely,  however,  as  e#  shall  see  later,  the  limitation  which 
sust  be  lnpoced  on  V  le  very  severe. 


he  shall  refer  to  (2.1.7)  as  Option  n. 

he  Dos  seek  a  fonmla  which  is  as  accurate  as  Option  II  and  which  Is,  at 

.  • 

the  sane  tins,  stable  for  all  r.  With  this  al*  we  expand  u(x,t)  In  a  Taylor 
aarias  about  tLe  point  (x^,t^)  obtaining 

\  ■  ”0  • *  bS»  *  |  t’v,  ♦  •  •  • 


•  5kSt  - 1  k\« 


Using  (2.1.2a)  and  (1.1)  we  have 

u,  ♦  hHu  ♦  ♦  $h2fc 


(2.1.11) 


1  ♦  hH  ♦ 

2r 


where  the  leading  tens  of  the  difference  between  the  right  and  left  cambers  of 

(2.1.11)  la 


(2.1.1 2)  ~  K^t  -^'xxx  . 


1  ♦  hH  ♦  ± 

2r 

evidently  tin  leading  tern  of  the  local  error  for  (2.1.11),  which  wo  stall  call 

2  t 

Option  III,  is  0(h  k  ♦  ha).  Uoreover,  since  all  coefficients  in  (2.1.11)  ara 
positive,  the  foraula  Is  stable  for  all  r. 

In  sons  cases  it  is  convenient  to  replace  by  f^  In  (2.1.11),  obtaining 


(2.1.13) 


^  ♦  hHu  ♦  ^  «4  ♦  -£b2f4 

; - ' 


1  ♦  hH  ♦ 

2r 


where  the  numerator  in  (2.1.12)  ie  changed  by  a  tarn  whoaa  abaoluta  valua  la 
§h*fcf^  ,  eta  re  f ^  *  f  t  ♦  f#i^.  Sines  this  tern  la  0(h*k),  the  order  cf  tha 
local  accuracy  of  (2.1.13)  la  the  sane  aa  that  of  (2.1.11).  be  will  refer  to 
(2.1.:  3)  aa  Option  III'. 


f 
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x-A~fa 


z*A  x**<b 


t-(j*l)k 


t-jk 


__ 

» 

.-JL.  I- 


Typical  Boundary  Canfii^patlon  -  Right  hand  boundary 

>lfura  S 

tt«  new  consider  the  treatment  of  (2.1.Zb).  lor  the  eir^l**  representation 
*e  replace  u^  by  (u^-  u,j)h  ^  obtaining 

(2.1.14)  ^  «/»  ttj  . 

It  is  easily  seen  that  the  leading  tani  of  the  difference  betaken  the  r\&it  and 
left  meber:  of  (2.1*14)  la 

-£  u 
2  xx 

«o  that  (2.1.14)  has  the  ea»e  accuracy  a«  Option  I.  Actually,  of  couree,  (2.1.14) 
le  a  special  case  of  (2.1.5)  wi*re  H  ■  0. 

lornula  (2.1.14)  represent#  the  only  xetLod  ahlch  we  uaed  for  treating 
(2.1.2b).  Psrhape  a  better  procedure  would  Un  been  to  eonaider  (IqJq)  ** 
an  Interior  point  and  to  uae 


(2.1.15) 


*1  ^  ®5  • 


Here,  the  leading  tern  of  the  error  is 


i 


Hi 


»*y  f  ■  0g!B  W  f\  >ei>  .  leg  +*m 


„  — r 

'^k.-  if*  V^t4  ijfii  *•  -3, 


•x-f  s 

■iikU 


t,  1,1  ch  It  0(h5). 

If  m  do..  not  td*  tt  ttt  tit  point  «“*»  “  <***  "* 

„  ^  .  tmrnl*  «  U.  u-  orttr  ol  ooconc,  t,  utUl  tto  I™**  <*«> *»> 

A.  formula  It  otol«d  f-  l*-MU  >»  1>UU»  *  ' 

ts  c°t 


(2.1.10) 


»,  ♦  £  *.  *  j*o 

"b 


A.  l.«lln c  ttnt  of  tb.  .rro,  1.  Oth5).  A.  -for.  th.  .cum*  .»ld  ».  JTC 
carved  If  «•  replaced  fQ  by  f4  In  (2.1.16),  obtalnin* 

r^*Zp^*  2h2f4  # 


(2.1.17) 


"b 


1  ♦ 

2r 


IS 


2.2  forward  difference  cat  hod.  Tor  the  derivation  of  this,  and  otter 
difference  aquations  it  will  be  convenient  to  consider  a  typical  configuration 
of  points  such  as  shown  in  figure  It., 


t-o*i)k 

.  ,  I 


i  I 


x*ih 


x-(i*l)h 


Typical  Configuration  -  Interior  Points 
figure  h 


Va  assume  hare  that  1  *i*M-l  ,  but  that  j  nay  equal  sero  in  which  case  the 
point  (x^  t^)  would  lie  outside  of  the  region. 

Probably  the  simplest  difference  equation  is  obtained  by  replacing  by 
k-*  and  by  replacing  u^  by  (ujSUj-2^)  h-^.  These  are  familiar  finite 
difference  representations  of  the  respective  derivatives.  However,  for  this 
wethod  and  tar  the  o there  given  below  we  shell  present  a  derivation  based  on 
wethods  of  numerical  quadrature.  Although  the  actual  finite  difference  form las 
com  out  the  same,  neverthelaaa  it  seewa  einpler  to  obtain  the  error  estiMtss 
by  the  wet hod  which  ee  use. 

3r  (1.1)  and  the  rectangle  rule  for  integration  we  have 


share  f  *  •  f  ♦  f  u, 
t  t  »l 


ktu^f )  ♦  V<tt«t*ft>» 


(2.2.1) 


/.:•>  -j: 


t„vk 

0  (u_  ♦  f) dt 


Using  Taylor's  seriss  we  also  obtain 
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(2.2.2)  b^  »  («*♦  Bj  -  ^  *\xxx  *  **#  # 

Substituting  in  tbs  previews  squstioo  we  gst 

(2.2.3)  Bj  s'  (l-Cr)Bp  ♦  r(a^  Bj)  ♦  kfQ 

wterw  ths  lsading  ten  of  the  difference  between  the  right  end  left  ■eeibers 
of  (2.2.3)  is 

(2. 2. It)  k  |(g  “  i^’Sotxx }  * 

Since  the  number  of  tine  steps  needed  to  reach  a  given  value  of  t  is  propor¬ 
tional  to  k~^  it  seeae  reasonable  to  tabs 

(2.2.5)  (  2  “  12  )'*xxxx 

as  a  Measure  of  the  local  error.  It  is  well  known,  (See  for  instance  [  18]) that 
k  oust  satisQr  the  condition 

(2.2.6)  k*Jh2 

in  order  for  the  metrical  procedure  to  be  stable.  Thus  the  local  error  for 
the  Forward  Difference  method  is  0(b  ). 

Of  course,  such  an  analysis  doss  not,  by  any  weans,  indicate  that  the 
accuracy  in  the  large  is  of  second  order.  She  question  of  the  order  of  con¬ 
vergence  in  the  large  has  been  studied  for  s  special  class  of  problems  ty  Juncosa 
and  Toung  [l2].  For  an  analysis  of  the  situation  in  the  large  for  ths  case  where 
u(x,t)  is  assumed  to  possess  certain  partial  derivatives  in  IMS,  see  Milne  [l7]. 

For  B  T  problems  with  t,$,f  vanishing,  ths  convergence  of  ths  sequence  of 
solutions  of  ths  Forward  Difference  method  to  the  exact  solution  of  ths  differen¬ 
tial  equation  has  been  proved  by  Lsutert  for  rf  l/4  ,  under  the  assumption 
that  g(x)  ie  piecewise  continuous,  that  is,  continuous  except  for  a  finite  number 
of  finite  Juapo,  and  that  g(x)  has  a  one  sided  first  derivative  everywhere. 
Hildebrand [9]  proved  convergence  for  0 <rg}  under  the  assumption  that  g(x) 


t 


Is  piecewise  continuous  and  baa  bounded  variation.  Mora  Mnrt  restriction* 
were  imposed  on  g(x)  far  the  earns  r*  |  •  Juncosa  and  Tenon*  [l_l]  proved  con- 
rsrgence  for  r  <$  under  the  assumption  that  g(r)  is  piecewise  cootinoous. 

The  numerical  computational  procedure  is  extremely  elmple.  Since  one 
knows  u^q,  0  •  i  ,  one  can  compute  u^A,  uaiag,(2«f.3)  which' 

for  a  general  point  becomes 

(2.2.7)  -  U-aO’Hj  ♦  r(^hj  ♦  VlJ*  ♦  kfA  . 

for  B  V  problems  Uq  a  and  u^  jars  glean  and  we  can  proceed  at  once  to  the  next 
roe.  For  N  D  problems  ee  can  compute  Uq  ^  directly  from  either  (2.1. 5^2.1. 7) 
or  (2.1.13)  and  ee  can  uae  (2.1.1b),  (2.1.1$)  or  (2.1.17)  to  find  w^  Having 
thus  determined  all  values  of  w^  ^  we  proceed  to  compute  values  of  then 

■ere  it  not  for  the  condition  (2.2.6)  on  the  sise  of  k,  the  Forward 
Difference  method  would  be  ideal  for  use  on  computing  machines.  Unfortunately, 
when  one  attempts  to  improve  the  accuracy  ty  reducing  h,  aay  by  a  factor  of  2, 
then  the  number  of  time  stepe  is  increased  by  a  factor  of  U,  making  the  total 
work  increase  by  a  factor  of  S.  More  generally,  one  can  assert  that  the  work 
Increases  with  the  third  power  of  b~l,  or  equivalently,  as  the  third  power  of  M. 


large  l8  0(h2)  proelded  k  =  0&  \  log  hj).It  i«  not  knore,  bo~ewr,  A.tlwr 
thlfl  result  can  be  extended  to  Include  cases  where  k  »  0(b). 


It  is  eoaetiues  convenient  to  replace  £j  bjr^  **  (?*3*2**  ^  addltloMl 

error  which  1.  Introduced  1.  proportional  to  k*  and  hence  tb.  expressloo  In 
brace,  occurring  In  (2.3-3)  is  Increased  by  a  ter.  proportional  to  k.  Thu. 

If  -  tends  to  tero  like  h,  the  expression  in  the  brwe.  would  be  0(b)  instead 


of  0(h2).  , 

xoother  .edification  Of  (2.3.2)  1«  obUlned  bf  replacing  the  Urn  ijf0*  tz\ 

by  •K*q**2)  }*  ThU  cb4n«*  ■“*  **  pr°<r“  pr*p*r*d 

for  the  Ordeac.  For  .UplUity.  we  as«u~  here  that  fU.t^ia  •  ***“«  of 


1 


i 
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a  ilon* .  Expanding  ln  a  Taylor  aariaa  about  tha  point  (x,t  ),  where  t^t^at^) 

ee  have  »  a  ♦  ^  ku^  ♦  g  k*t»tt  ♦  •  •  • 

"b  -  n  -  }  k^  ♦  J  k\t  ♦  .  .  . 

fa  '  7  *  *  “t  *  8  k W  *  *  * 

f0-7-*kft48k2ftt"*-* 

— *  «.  12^ 

■hara  f  -  f(u),u  *  u(x,t)  .  It  is  evident  that  J(fg*  fg)  *  *  4  g  k  ftt  4  •  •  • 

and  i(^,a  t^)  -  «  4  r  k  ntt  4  *  *  * 

Therefore  f  (H1^*  Sg)  )  *  7  *  0  k^tt^  4  *  *  * 

«d  Kf0«z)  -  «J(  V»2> )  *  i  k'f't  -  J  k2*^ t  •  •  •  •  •  «<**>• 

fy  (2.3*2)  and  (2.3.3)  it  is  evident  that  tha  order  of  accuracy  has  not 
been  changed. 

For  a  general  point  (2.3.2)  aay  be  written 

(2.3.U)  *if j.!  ■  «ltJ  *  *  Vi,j  *  "i.i.j.1  *Vi,iJ  _ 

Two  alternative  numerical  procedures  are  available  for  using  the  Crank-^**^  ' 
licolson  nethod.  The  first  is  an  iterative  aethod  wherein,  given  all  values 
of  uj  j,  OaitM,  one  chooses  a  fir.it  approximation  0*1  and  pro¬ 

ceeds  to  iterate  using  the  formula 

(l.J.S)  us; j.l  •  “jljf  *lfJ  *  2(1*0  K.1.J  *  *  a  1-1. J*1 

*  iufcykj  *  fni]]'  <->■!?; 

Here  the  parameter  “  known  as  the  relaxation  factor,  is  chosen  in  the  rang* 

1  4  »  <  Z  to  accelerate  the  convergence,  it  the  end  of  each  iteration  one 

obtains  (2* 1.5), (2.1* 7),  (2.1.U),  or  (2.1.13)  for  ^^uslng 

the  latest  values  of  j^ehen  appropriate.  Similarly  u^  ^  is  obtained. 

The  iterative  procees  it  repeated  until 

Hex  I  e(n)  -  m(n"l)  \  <  s 
Otttl  i#J*l 

there  «  is  e  preesslgaed  tolerance.  Vma  the  iterative  process  has  converged. 


(2.3.6) 

there  • 

one  proceeds  to  dstendm  the  values  la  the  next  roe  in  the  seme  way. 


The  choice  of  the  relaxation  factory, which  will  /laid  the  futMt 
can  easily  be  cospoUd  for  linear  BV  problems •  It  is  shown  in that 

wharw  I4  is  the  eigenvalue  of  the  (M-l)x(ll-l)  *trlx 

/opoo  ...  0\ 

/  po  f  0  ...  01 
I  0  p  0  p  ...  0  ) 


\0  0  0  0  ..p  0/ 

where  p  “  It  is  easy  to  show  that 

(2.3.8)  P*^bco,g* 


(2.3.8)  *  "  C01  M  ‘ 

Let  us  now  assume  that  k  =  0(h)  and  that,  as  a  result,  r  is  propor¬ 
tional  to  M.  Since  for  large  II  we  haws  ,  we  can  show,  as  in  [&], 

section  4,  that  the  rate  of  convergence  of  (2.3.5)  varies  as  lf*j  hence  the 

number  of  iterations  increases  as  II*.  Since  the  number  of  points  oo  each  row 

3/2 

increases  with  K,  it  follows,  that  tha  work  per  row  is  proportional  to  M  , 
Uoreover,  since  the  number  of  time  steps  increases  as  K,  (instead  of  M2  as 
with  the  Forward  Difference  wethod)  we  conclude  that  the  total  time  increases 
as  M5^2with  the  CrarJc-flicolSon  method  compared  with  for  tha  fcrward  Differ¬ 
ence  method.  Thu;  the  Crank-Nlcolson  method  waild  appear  to  be  miperlor  to 
the  Forward  Difference  method,  et  least  for  large  H. 

The  advantage  of  the  Crank-Nlcolson  method  can  be  increased  still  further 
by  using  s  non-iterative  procedure.  The  effectiveness  of  this  procedure  de¬ 
pends  on  the  fact  that  for  a  linear  system  involving  a  Jacobi  matrix,  l.e.,a 
matrix  whose  only  non -ter  o  elements  occupy  the  main  diagonal  and  the  two  diagonals 
adjacent  to  the  mala  diagonal,  the  solution  can  be  obtained  very  easily.  This 
fact  was  obsarvwd  by  l.H.  Themes  in  unpublished  notes  [2]]and  was  applied  to 
parabolic  equations  by  ftruct,  Peaceean,  Rachford  and  Rice  [2)  .  Tlais  following 
their  discussion,  let  us  suppose  we  have  the  system 


r 


(2.5.9) 


(  Vl  4  Cl*2  "  dl 

Aixi-1  4  BiXi  4  *  di 


£  iliM 

3*?  . 


ViVi '  V»  ‘S 

If  on*  carri**  out  th*  faniliar  Oauss  eliadaaticn  procedure,  elimi¬ 
nating  th*  element*  be  Lem  the  udn  diagonal  of  th*  matrix,  and  than  uaaa  th* 
method  of  back  substitution,  on*  obtains 


(V'i 

" 

(2.3.10a) 

1*1  •  ’I  -  Vi.1 

1  tfi  i  *  H-l 

ah*  ra 

( h  •  “A 

(2.3.10b) 

dl  ~  Vk-i 

Bi  “  Aibi-1 


and 


(2.3.10c) 


bo"° 


Bi  -  Aibi-1 


1  9  1  f  5-1. 


The  stability  of  tlu  confutation  procedure  has  been  investigated  by  T  items* 

in  [2]. 

No  rare  than  about  3  or  4  tints  as  much  cork  per  roe  would  appear  to  b* 
required  with  this  method  as  with  th*  Forward  Difference  method.  Thus,  if  w* 
aasune  k  ■  0(h),  t>*n  th*  aaount  of  cork  required  increaa**  as  M*  Instead  of 
M^2  with  th*  iterative  procedure  and  with  the  Forward  Difference  matted. 

Unfortunately,  when  (1.1)  ia  non-linear  it  la  not  alwmy*  eaay  to  obtain 
linear  equations  and  still  retain  the  desired  accuracy.  For,  if  w*  us* 
(2.3.2),  we  find  that  th*  unknown  vain*  of  *  la  involved,  in  fsnerwl, 
is  a  non-linear  manner  in  th*  evaluation  of  th*  function  f(z,t,«).  Cb  th* 


r 


*  ■-+-  -*».*  --v  — -  • 


r  V. 


*  J  « 

*  * 


-  -r1 


l 
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other  hand,  if  cr.e  replaces  u^  ,+i  by  j  ,  one  obtains  a  Ism  accurate  f  araiula, 
as  discussed  above.  To  get  around  this  difficulty  in  an  analogous  situation, 

Bruea,  Paacesan,  Rachford,  and  Rica  (2^  used  an  iterative  procedure  together 
with  the  ncn-iteratlva  aethod  for  a  airing  the  linear  equations.  Although  this 
procedure  sight  appear  to  defeat  the  original  purpose  of  avoiding  all  itera¬ 
tions,  it  is  reported  that  convergence  is  obtained  after  a  very  few  iterations. 

It  eight  not  be  too  unreasonable  to  praassi£i  a  fixed  number  of  iterations, 
independent  of  h  and  k,  to  be  carried  out  at  each  tisa  step.  Thus,  the  Ueun 
aethod  for  solving  ordinary  differential  equations  is  the  ssae  as  the  aodified 
Euler  netted  (see  fur  instance,  Uilne  [l?])  except  that  one  performs  only  one 
iteration.  The  local  error  for  the  Koun  aethod  is  sonewhat  larger  than,  but  j 

I 

of  the  saae  order  of  magnitude  as  the  aodified  Euler  aethod.  i 

Another  possible  solution  to  this  difficulty  would  be  to  use  a  difference 
equation  involving  more  points.  However,  care  aust  be  taken  to  svold  formulas 
which  night  lead  to  unstable  procedures,  and  also  to  avoid  obtaining  a  matrix 
which  was  not  a  Jacobi  matrix. 

j 

2.4  Pur  ort-r  ranks  1  method.  In  [5]  Dufort  and  franks 1  obtained  a  difference 
equation  which  ia  stable  for  all  values  of  r  and  which,  st  the  saae  tine,  yields 
explicit  formulas  for  values  of  u  in  the  new  row  in  tens  of  values  on  previous 
rows.  aay  derive  tie  for  aula  by  use  of  nuaerlcal  quadrature  as  folios. 


ft  *« 


Using  Taylor's  series,  as  before,  we  f*t 


f  2 

(2.4,2)  Wj  -  ^  ♦  2^^  ♦  *5  -  ♦  fj  ♦|-aki£,a»na 

4^<Vrt  ♦'!>}  •  • 

If  we  were  to  neglect  tbs  expression  in  the  braces  and  all  b  jher  order  tern 
in  (2.4.2)  we  would  obtain  Richardson's  method  [20] ,  which  Is  unstable  for  all 
values  of  r  and  is  thus  unsuited  for  numerical  calculations  unless  veiy  special 
precautions  are  taken  [is}.  Nevertheless,  at  this  point  we  remark  that  the 
local  error  for  Richardson's  method  is 

(8.4.J)  k[-£'w*  sr<wft>/  • 

In  order  to  obtain  a  stable  method,  Dui'ort  and  b'rankel  replace  Ug  by 
i(u4  ♦  «  )  in  (2.4.2)  and  obtain,  after  solving  for 

(2.4.4)  *2  •  \  4  ut  (U1  ^  U3}  4  I S?  f0 

♦  ife  £•  V  0xx« 4  s"  (\xt 4  v "  5^  ^  j 4  ** 

Now  if  k  •  O(h^),  it  is  evident  that  tha  expression  in  the  braces  is  0(h*).  On 
the  other  hand,  if  k  ■  o(h),  then  the  expression  inside  the  braces  need  not  even 
tend  to  s«ro  with  h.  Consequently,  although  stability  is  assured  for  ill  h  and  k 
for  convergence  the  ratio  l^h  must  tend  to  sero  with  h,  (see  [sj). 

Neglecting  remainder  terms  we  here  for  e  general  point 

(2.4.5)  ^  \yl  ♦  ife  (“l*l»J  4  “1-1, P  4  1<?  fl»J  * 

evidently,  to  compute  u  .  one  needs  the  values  of  u  on  the  >-th  and  (>»l)st 

1,3*1 

rows,  it  the  beginning  one  has  only  the  values  far  J  •  0.  To  obtain  tbs  values 
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of  u  for  J  •  1  another  procedure  must  bo  ased.  This  paint  it  discussed  further 
in  Section  4. 

The  convergence  of  the  Cd'ort-Frankel  atW  for  *  certain  clegs  of  prob¬ 
lems  lies  been  given  ty  Juncosa  £l4^a 


2.5.  Extrapolation  to  zero  grid  alee.  For  the  proceea  known  ee  extrapo¬ 
lation  to  zero  grid  site,  proposed  bj  L.  tlcherdson  ^20 J,  one  assumes  that 
the  error  ag(xyt)  »  n^(z,t)  -  u(x,t)  la  proparticnal  to  h*.  Here  u(x,t)  Is 
the  exact  solution  of  the  differential  equation  and  *^(x,t)  la  the  solution  of 
the  difference  equation  with  h  “  l/H.  Usually  certain  restrictions  are  assumed 
for  k.  Even  though  e^Ujt)  night  vanish  like  K  ,  as  la  indeed  true  for  certain 
cases  as  shown  in  [12] ,  [is],  nevertheless  it  does  not  necsssarlljr  follow  that 
.  (a.t)  is  exactly  proportional  to  M  .  This  can  be  easily  seen  by  considering 
tin  esse  where  e^(x,t)  "  if* sin  M. 

If  were  directly  proportional  to  if* ,  ana  could  determine  u(x,t) 

in  terns  of  u^fxjt)  and  u2^(x,t).  Indeed  we  would  have 

(2.5.1)  u(x,t)  -  ^(x.t)  *  ^  "  %(*»r)]  • 

If  one  solves  the  difference  equation  for  three  different  we  oh  sizes,  ob¬ 
taining  u^(x,y),  u^(x,y),  and  u^(*»7)  *•  em  perform  a  siapls  test  f dr  the 
assumption  that  e^(x,y)  is  proportional  to  K  .  Indeed,  if  the  assumption 
were  valid,  we  would  have 


(2.5.2) 


•w(*»y)  j 


By  cor put log  the  ratio,  p,  appearing  la  (t.U)  one  can  determine  the  proper 
extrapolation  to  sere  grid  else,  provided  sweh  exists,  end  also  estimate  the 


n 

digram  at  conference  of  the  sequence  of  solution*  of  tho  difference  equation. 
U,  oat  MCMM17  to  compute  a  itself  since  2a  eppeere  in  (2.S.1).  *e  wish 
to  eepheoino  aoce  00 re,  that  the  considerations  of  this  section  ero  net  rifar- 
oos. 
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s.  specific  g£fifeltaa  Sscald&ofl 

In  order  to  obtain  sosm  Indication  as  to  the  validity  of  the  non-rigorous 
oat  he  optical  considerations  of  the  previous  section,  we  havo  solved  sossi  sajqple 
problems  on  the  Ordvac  cocputer  at  Abordeen.  be  chose  problems  involving  either 
boundaiy  conditions  (2.1.1)  or  (2.1.2)  and  either  linear  and  non-linear  differ¬ 
ential  enactions.  Actually,  as  will  be  seen  in  the  next  section,  the  program 
which  was  prepared  for  Ordvac  is  capable  of  handling  entire  classes  of  problems 
rather  tl>an  no  roly  the  specific  problems  described  belcm. 

3.1.  Boundary  value  oroblea  -  linear  differential  equation  (BVL).  Here 
we  selected  a  problem  for  which  an  analytic  solution  is  readily  available.  Tim 
equations  are 

iu  «u  ,  0  <-x  <  l,  t  >  0, 

t  xx 

u(x,0)  «  g(x)  ,  0  <  x  <■  1 

u(0,t)  -  u(l,t)  •  0  ,  t*0. 

The  analytic  solution  for  tlas  case  g(x)  ■  1  is  given  by 

CO  ^  2 

(3.1.?)  u(x,t)  -  \  £  J  sin  a»  t  "n*"  t 

r.-l 
n  odd 

be  flail,  by  solving  the  difference  equations  with  various  mesh  sizes, 
attempt  to  determine  ths  rate  of  convergence  for  both  methods.  For  the  Crarvk- 
Nicolson  msthod  we  slall  let  k  -  0(h)  ss  h-eO,  evsn  though  the  theoretical 
results  have  been  proved  only  under  stronger  restrictions  on  k,  £isj,  in  an 
attaint  to  deUroins  whether  a  genersliiatioa  Is  to  be  expected. 


3.2.  Boundary  value  problem  -  ncn-lineec  MXgrgfftjM  taaslla  I S5CU 
Tbs  foliating  p  rob  lea  described  by  Hart  re#  [?3  appears  to  bs  typical  of  this 
class  of  pi  oblaaa. 

0*x.<l,  t>0 

(3.2.1)  1  u(x,0)  » 0  0  *  x  <  1 

u(0,t)  •  u(l,t)  *0  t  feO 

Here  u(x,t)  represents  tita  temperature  in  a  substance  where  tiers  la  generation 
of  heat  increasing  exponentially  with  the  teaperature.  *hen  (  exceeds  a  certain 
critical  value  tiers  is  no  steady  state  solution  and  the  solution  increases  be¬ 
yond  all  bounds,  be  shall  not  study  this  aspect  of  the  problem  but,  rather, 
will  use  a  fixed  value  of  0.8  for  f ,  which  is  below  the  critical  ralna.  We 
shall  solve  the  problem  by  the  larward  Difference  and  Crank-Hicoleao  method 
for  several  values  of  h.  Solutions  for  various  values  of  y  have  been  obtained 
by  Kart  rue  using  a  differential  analyzer  and  tha  results  agreed  -*uite  wall  with 
experiments. 

3.3.  Norm  I  derivative  problems  (~M)L  end  NDH) .  Here  »e  consider  b<mcx'iu'7 
conditions  (2.1.2)  with  both  the  linear  equation 

(3.3.1)  u  -u  0  <  x  <  A,  t  >  0 

T*  XX 

and  the  non-linear  enuatioa 

(3.3.2)  \  ♦  •‘Vtt  ,  0  <  x  <  A,  t  >0. 

The  initial  conditions  are 

(3.3.3)  f (x)  -  e  0  III  A 

where  in  tie  caeee  which  we  stall  eons  ids r  c  •  .051.  Actually,  in  our  eeejut*- 
tiona,  *e  stall  take  for  our  initial  functioo  not  g(x)  ■  .031  but  a  eat  rf 
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values  obtained  for  a  later  ralue  of  t  using  the  analytic  solution  of  the 
linear  ,-roMen  (rue  for  instance  Ms  page  22) .  The  values  are  given  is 
TaVle  I.  Even  in  the  mm-llnear  case,  these  valves  arc  sufficiently  accurate 
since  the  non  linear  tern  ie  negligible  for  eju.11  values  cf  t.  This  adjust- 
nent  in  the  initial  conditions  was  cade  to  eliminate  ary  pot  viola  effoct  of 
the  singularity  caused  1 13  the  discontinuity  of  u^  r<*an  t  *  0.  It  io  by  no 
means  certain  tl*at  such  a  precaution  would  te  rxceesary,  and  for  the  linear 
case  at  least,  it  should  be  possible  to  dc  Lerralna  analytically  the  effect  of 
tie  singularity. 

Problems  of  tids  type  involving  (3.2.1)  are  designated  ®L  and  those  in¬ 
volving  (3.3.2)  are  dusign.it.ei  "*>!.  It  or  .3.'  probloas  tin  solution  u(x,t) 

re.-re;  unta  the  te:  perulure  in  a  propellent  uliere  the  reaction  rate  ia  deter- 

-l/u 

njiiud  ly  tie  Arrhenius  ter#  e  .  The  prcpellant  is  assnnod  to  be  a  sard- 
infinite  slab,  0  *  x  <00.  Tie  propellant  ia  leaUd  try  a  gas  tt  tearentur* 
u  i.l. ich  occupies  the  region  x  <■  0  and  whose  influe. aco  on  tha  propellant  de- 

ia 

iKT.ds  u«»  the  1  jai  transfer  coefficient  h.  In  order  to  be  able  to  use  finite 

difference  methods  »e  consider  a  slab  of  finite  thickness  A  and  replace  the 

condition  Lin  u  •  0  by  the  ctnditiun  (2.1.2b).  rurther  details  on  the 
x-00 

vs  Seal  rotivstian  for  these  protista  c^y  be  found  in  [8 J. 

In  outer  to  represent  *ore  accurately  tiie  ptyeical  situation,  it  ms 
conutir.js  necessary  to  introduce  "shutoff",-  tiat  le  after  a  certain  value  of 
t,  is  greatly  reduced.  However,  te  :lall  not  do  this  in  ear  study* 

be  shall  also  not  consider  tl«e  (letsrrirwtion  of  the  "ignition  tine"  which 
ia  the  tlra  required  for  tta  value  of  a(x«t)  to  incroaso  above  a  preaaaignad 
value. 
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3U11  S3SUIW  the  follow  1a:  *»!«•■  for  thc  P*”*"*^” 

4  -  75,000 
e  •  .051 

«  •  .18 

f  -« 

H  -  2*10 

Prior  to  1952  a  number  of  calculations  were  dene  oo  this  prcttsn  al  the 
Abordcon  Proving  dround  using  the  Bell  Relay  Computer  and  the.  Bniae.  The  For¬ 
ward  Difference  ret  hod  was  used  with  r  -  £  and  several  values  of  h.  Is  shall 
aake  use  of  core  of  tlese  results  in  Section  5. 
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4.  Description  o£_  £hp.  Prcxran  frr  bHit 

To  run  a  particular  problea  with  ti/e  pro^r*n  as  precartly  coded.  It  Is 
necessary  to  specify  three  major  conditions)  l.e.  (1)  tie  dlffererte  »*,u.itlce» 
to  b*  used,  (2)  tLe  bcund&ry  conditions,  and  (3)  the  r.eceiJury  paransters  such 
as  h,  r,  etc. 

All  of  tie  difference  situations  and  boundary  conditions  described  in  Sec¬ 
tion  ?  liave  been  ceded  ujL  any  problem  involving  these  can  be  computed  using 
tin  program.  Sxcept  for  te  CKT  nethod,  tho  interior  pcint*  of  the  (j^ljtt 
roe  are  computed  followed  by  tl«e  coeputatiun  or  insertion  of  tie  boundary  values 
of  tint  roe.  Ti.1*  perrdtc  the  boundary  conditions  to  lie  a. icily  codified.  The 
*p.endix  contains  the  detailed  enuatione  used  in  the  present  coding. 

Unleee  othcreiea  instructed,  the  r  *. chi  re  .ill  auturae  tl-.t  tl.e  differential 
tr.|Vr.tian  lc  li.T&ar.  In  tl.e  rwr-lircar  case  it  is  necessity  to  instruct  the 
-aci.ino  to  caapvl*  the  r on- linear  tern. 

lor  tl<c  flia  eethed,  since  the  j-th  arid  (>*l)st  r«u  i  r«r  necesstry  to  coo* 
pu'.e  tl.e  (j*l)ct  rt*.  (;c«  appendix),  it  vac  r^citrary  to  ccrpute  tlio  first  roe 
ly  a  different  »thod.  »«  used  tie  S3  ecth'xi.  Tbs  rueher  of  rwe  necessary 
to  le  C deputed  ly  the  rD  fttthod  defends  upon  r,  ..lush  involves  k,  tie  i.rcre- 
"ciit  in  tieo.  Tids  is  noted  in  the  A^iendix. 

lor  tie  CM  sothed,  HUH  case,  men  tie  valors  of  tl.u  function  becaee  larx**, 
the  iterative  echeee  would  not  converge.  A  future  was  added  ui  the  prograM 
‘.beieb7  If  tie  rusher  of  iteration*  eiceeded  a  .liven  for  the  (j*l)et  rew, 

Mu  coding  wce.ld  revert  buck  to  the  tO  retl.ai  and  co.-jn.ince  c>sj,A.t*.tiun  f roa 
tie  >*th  roe  urine  r\  'here  r  •  r'«2®  and  0  <  r*  <  J. 

SI  ait  off,  another  future,  can  be  describee  a t  reducing  tl.e  gas  terpirature 
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<V  or  shutting  off  tbs  gu  flow  after  tin  tg  (eee  Section  3).  after  com¬ 
puting  row*  (TO  only),  the  value  of  can  be  clanged  for  the  ( Js4l)*t  roe 
and  the  retaining  row*. 

The  pars  raters  that  aust  l<  choeen  are  I 
b  -  ta/A 
r  -  k/h* 

S  switch  -  which  is  the  stopping  procedure  where  3  •  stops  the 
enchine  after  a  roes  and  S  ■  S,,  stops  macliino  elan  mtx  |u.  ,  | 

(C  *#J 

exceeds  or  falls  below  a  given  tolerance. 

Besides  these  paraaeiers,  certain  specific  cases  will  require  a  specifica¬ 
tion  of  other  inf  creation,  such  asi 

J^(CN)  -  as  explained  above 

«  (Cb)  -  the  relaxation  factor.  If  u  is  not  specified,  *>  *  1 
li  assoasd 

y  (BVN) <  (.8  was  used) 

H  (TO)  $  (2*lo”®  was  used) 

Bg  (TO) »  (.18  was  used) 

A  (TO)  <  (73,G00  was  used) 

j  ,  u  ^  (TO)  -  (shutoff) )  replaces  u  by  u^  after  coaputing 
a  (  A  8 

rcwe. 

The  non-linear  tent  was  cccjxited  by  a  floatiry;  point  sub  routine  with  the 
exponential  factor  being  cceputod  by  taking  the  first  15  terws  of  tha  series 
expansion  of  ez.  Initially,  the  exponent  was  nonaalitedj  i.e.  placed  in  the 
fora  a*2*  where  J  t  |a|  1,  or  e  •  0  where  appropriate,  (tor  the  BV  case*, 

«  -  4  since  ell  values  were  scaled  by  2~*.  tor  TO  cases,  the  preaenoe  of  the 
exponent  (-  })  necessitated 


-IP***.  p**-. 


so 


-4  * 

the  nonnall*in£  prtcvdtre.)  2  •  «a*  co  puted  foliated  by  s  narunlixetiaa 

of  e*.  Sine*  mu  desired,  the  scaling  problem  »as  solved  by  the  following! 

.'(*)  -  -  <e*)2* 

Therefore,  it  was  neereeuy  to  Stuart  a*  a  tires  in  cob putlog  the  ncn-I  Irear 
term. 


5.  Antlysls  of  Bssulta 


In  Tables  II  -  IX  we  give  the  values  of  u(x,t)  for  the  respective  prob¬ 
lems  as  obtained  by  the  various  finite  difference  methods.  Sealing  factors 
are  indicated  an  the  tables.  In  each  case  the  columns  headed  *D($)  and  *D(J) 
contain  results  obtained  using  the  Forward  Difference  method  with  r  •  |  and 
r  “  1 /4,  respectively.  Columns  headed  CM(I)  and  CMt(l)  contain  the  re  suite 
obtained  by  using  the  Crank-Nicolaon  method  using  the  iterative  and  noo-iter- 
atlve  numerical  procedures,  respectively.  Kor  different  values  cf  li  ■  k/ b, 
the  following  values  of  r  are  used 

JL  -X_ 

8  1/2 

16  1 

32  2 

64  4 

The  column  headed  CNH(II)  contains  results  obtained  using  tie  Crank-Mlcolson 
method,  non-iterative  and  with  tie  following  values  of  r 

AZJ1  -JL- 


8  1 

16  2 

32  4 

64  8 


The  value  of  t  can  U»  obtained  from  the  roe  number  by  multiplying  the  roe 
minbor  ty  Thus  th®  ro"  n’J8**r  t/kg  where  kg  is  the  vslue  of  ' 

k  rorresporvllng  to  U  •  8  and  r  ■  1/2.  The  column  headed  *Kniac"  gives  results 
obtained  on  the  Snlac  using  tie  Forward  Difference  method  and  r  *  1/2. 

Corresponding  to  each  of  the  rem  numbers  4,  8,  12,  16,  and  20  and  to  the 
indicated  values  of  x/A  »*  give  tbs  following  results,  when  available, 
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Hera  u^  represent*  the  solution  of  the  difference  equation  being  considered 

with  h  -  A/U.  Also,  *16  -  o16  -  Ug,  A  -  B  “  #S2  "  ~  *16»  C  '  *&/*&*> 

and  D  -  e  -  u  -  u,„.  In  tie  brackets  »e  give  the  extrapolated  value.  The 
64  64  oZ 

estlowited  value  of  p  used  in  the  extrapolation  is  given  at  the  bottoa  of  the 
colunsi,  (seo  section  2.5).  for  tlxs  extrapolation  e^  is  used  alien  available, 
otherwise  eJ2  is  used.  An  asterisk  in  place  of  e^  indicates  that  the 
c Deputed  value  is  not  regarded  as  slgnlf leant,  usually  because  e^  is  too  snail. 

An  asterisk  near  an  extrapolated  value  wa-'.s  that  the  value  obtained  from  tlie 
use  of  the  scaliest  eesh  size  is  jiven  in  place  of  the  extrapolated  value. 

5.1.  poundary  value  linear  (P7L).  In  this  care  the  results  obtained  fro* 
tl*i  Crank-Nicolson  s*tnod  »ith  the  iterative  and  tlie  non-itaratlve  mxa»rical 
procedures  were  identical.  Since  the  rstlos  %  /' •jyj  to  4  it  wbuld 

s 

appear  tiat  eg  -  odT2)  and  tUt  extrspclaticn  to  sero  grid  size  based,  on 
a  •  2  might  give  goad  results.  Indeed,  if  one  uses  this  procedure  one  obtains 
altost  identical  results  for  all  *»thods  and  very  close  agree nent  with  the 
exact  value  aa  calculated  fro*  (3.3.2).  U  note  thst  the  values  obtained  fro* 
the  forward  Difference  set  hod  sith  r  -  5  and  2  -  52  are  considerably  different 


f  roc  the  exact  values  end  that  the  extrapolation  process  has  greatly  lap  roved 
the  accuracy.  '«*  also  remark  that  tor  this  particular  problem  the  Crank- 
Hlcolson  method  is  such  more  accurate  than  the  iorward  Difference  method 
although  the  extrapolated  values  are  practically  the  same. 

5.2.  Boundary  value  non-linear  (BVN).  As  in  the  linear  case  the  values 
of  are,  with  tha  kD  and  CS  methods,  close  to  4.  Extrapolation  to  sero 

grid  site  using  p  •  4  and  a  •  2  yields  results  which  are  nearly  the  eaaa  for 
both  tlie  ID  method  (with  r  •  J  «  sell  as  with  r  *  1/4)  and  for  CN(I).  On  the 
other  iiand  with  the  non-iterative  procedures  CVH(I)  end  CNN(II),  (where  the  • 
less  accurate  difference  equation  obtained  by  replacing  by  fQ  in  (2.2.2) 
was  used),  tha  values  of  e^  /e^  were  usually  between  5  and  5.5.  Thus  as 
ex;*ected,  the  order  of  convergence  was  less  than  for  the  other  methods.  Nev¬ 
ertheless  with  an  extrapolation  using  p  •  S  »e  obtained  fairly  accurate  results 
especially  for  C.1N(l). 

'••e  believe  tint  caabiiung  C.VN  with  a  few  iterations,  not  core  than  five, 
would  increase  tie  order  of  convergence. 

The  »*rro r  due  to  using  f  instead  of  f^  in  (2.5.2)  varies  approximately 
as  M  *  as  can  be  seen  froa  cosparing  rusults  of  CR(l)  end  CflK(I).  Thus  for 
tie  10th  ran  and  for  x  •  J  we  tiavs 


11 

CK(I) 

ChK(l) 

aKi)-cw(i) 

8 

72554 

72400 

154 

16 

72668 

72821 

67 

52 

72977 

72945 

54 

(ill  velmas  have  been  aultlplimd  by  10®.) 
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Explaining  the  differences  134,  67,  34  w  ot se rve  that  they  are  approximately 
in  the  ratios  4i2tl.  Hence  the  differences  appear  to  tend  to  aero  like  If”*. 

This  tends  to  confira  our  prediction  that  the  error  introduced  by  using  the 
less  accurate  difference  equation  is  of  order  M  \ 

5. 3.  ygnsaj  d_c^v.vtivp  Ucaar  ian.  In  both  tie  linear  and  the  non¬ 
linear  cases  we  are  ccncerred  with  tlirec  options  for  representing  the  condition 
(2.1.2a).  he  first  consider  the  result*  obtained  for  the  linear  case  for  Options 
I,  II,  III  ac  given  in  Tables  IV,  V,  VI  respectively,  as  in  the  boundary  value 
linear  case  it  did  not  matter  whether  one  used  the  iterative  or  the  non-iter¬ 
ative  procedure  for  tlw  Crenk-lficolscn  oathod. 

for  Opticr.  I  the  values  obtained  by  the  borwuri  Difference  method  for 
r  •  1/4  and  r  -  1/2  agree  closely  with  ea'  h  ether  and  with  those  obtained 
with  CNN(I)  and  CCi(II).  Based  oh  the  observed  value  cf  p,  extrapolated 
values  were  cbcair^J  as  indicated.  Of  course,  in  this  case  we  cun  conclude 
little  from  the  fact  t!ut  tho  extrapolated  values  obtained  from  the  different 
not  hois  agree  with  t.'h  ether,  iiisvover,  the  extrapolated  values  do  agree 
rather  closely  -lth  triune  obtaint-d  frea  results  using  Option  II  and  III,  and 
tho  latter  Uo  agree  with  e.«ch  other.  Wcieover,  the  difference  -  u^ 
obtained  lor  CMJ(l)  with  Option  III  was  very  scull,  and  this  indicates  that 
UG<  *5  probably  very  close  to  Uie  tn*e  value  of  u.  Thus  it  would  appear  that 
in  this  case  (2.1.5)  (Option  I)  is  a  reasonably  satisfactory  representation 
of  (2.1.2a)  provided  cne  uses  the  proper  extrapolation  procedure* 

finally  we  compare  our  result#  *lth  the  "exact*  values  obtained  under  the 
aarx^tlcn  that  A  is  infinite,  using  the  procedure  given  in  (6).  These  values 
are  given  on  T«ble  1/  in  braces.  Bear  the  line  x  •  A  our  values  will  not  agree 
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«u>  «i-«.  “  «“  ““  *  u  rtnl“-  B"'rth‘1'“  “* 

agr^uuit  1*  raUier  cloee.  . 

ie  ^  m  table  V  the  values  obtain*  on  the  ini*  using  Optim  U  and  . 

the  H>  nsethod  «ith  r  -  J.  Although  the  Enlac  values  and  the  Ordvac  values 
listed  under  »D,  again  «ith  r  •  agree  exactly  for  M  -  8,  exact  a«r«.«nt 
cannot  be  exited  for  other  value,  of  H.  The  reason  for  thi.  Is  Uat  A-64,400 
„9  u=ed  for  tb.  Sru.ac  results  Instead  of  A-73,000  as  used  «  the  Ordv.c  to 
the  Snisc,  instead  of  (2.1.14),  the  condition 
*  “Bel 

...  „.d  to  represent  (2.1.2b).  "«■»,  >*«•  *  l'»  *■  “*;t 

i,;roeuwnt  in  the  case  M  *  3. 

Near  the  left  hand  boundary,  p  apixars  to  be  rather  close  to  4.0  so  that 
w  use  extrapolation  based  on  a  •  2.  We  note  that  is  snail,  and 

this  indicates  that  uJ2  is  already  HUite  accurate,  lor  xj k  ”  $  the  use  of 
p  -  J  indicated,  and  for  x/A  -  1  ~  »•*  p  -  2.  It  appears  tUt  the 

largest  errors  occurred  at  the  right  hand  boundary. 

Sirjl.r  remits  ver e  o'  talned  Sci  Option  III  a*  located  in  Table  V. 

Of  course  for  Option  III  it  .as  possible  to  use  CNH(I)  and  thus  use  nany  fe-er 
titt)  .tep.  U»n  .ere  necscary  for  Us  Ior>ard  Difference  nnthod. 

It  sec  i*.sc»»Mo  to  suppose  licit,  ty  using  either  (2.1.15)  or  (2.1.17) 
to  represent  ror*  accurc’oly  the  right  hand  boundary  condition,  and  by  using 
Option  III  or  III',  the  uUo  »  uculd  be  nearly  4  ever^b^re  and  that  extrapc- 
lavicn  Used  on  e  -  2  could  Uve  been  used  for  all  points. 

5.4.  awl  WYhUa  H£n=^  iHO-  **  r,3UlU  eU#u,#4 


using  the  ID  method  lor  Option  II  indicate  that  extrapolation  to  stro  grid 
siM  using  o*2  for  x/k  *  0,  l/k  and  a  •  1  for  x/k  •  1/2,  1  should  bo  quit# 
accurst#.  Tie  saae  also  appears  to  be  true  for  our  Ordvac  results.  Of  course 
the  axtrapclated  values  obtained  in  th*  respective  cases  do  not  ajyee  because 
of  the  difference  in  t ne  boundary  conditions  at  x/k  •  1,  aa  oantioned  abovo. 
Tba  extrapolated  values  cbtairod  cn  the  Ordvac  using  Option  II  and  the  M) 
wethod  also  agree  very  closely  with  tiiose  extrapolated  values  obtained  using 
CH(1)  with  Option  III  It  is  difficult  to  sake  definite  statenents  as  to  the 
behavior  of  the  CSN  anthcd  except  tint  the  errors  sre  snail  aril  irregular, 
here,  of  course,  Option  III'  was  used. 

for  Opticr.  I  ti«  r;ilios  p  are  all  a-proxin  tely  2  with  all  difference 
e^uaUr-ns  as  in  t*«  linear  case.  Tha  errors  are  ecwwewhat  larger  and  the 
cxtr-ioolatud  values  acw«.'«it  less  accurate  t!an  with  Option  III, 

5.5.  :lich, n«  .tijrrs..  figures  5  and  C  rhew  the  approximate  tines  required 

to  c  depute  t  •  lOg  revs  by  the  vuricus  reticle  m  Loth  the  lihear  and  in  tie 

n  on- 1 1  near  ‘-.is-.-s  vt  >ly.  Tie  t..v>’  for  BV  prooleus  and  for  HD  problina 

*.-:e  r«arly  U*  an/..  Tie  ri a  for  tin  curves  c^rrecpcndlng  to  tiie  fD,  CH 

«.-.d  CHI  rnti.jdr ,  after  proper  correction?  have  been  aaJe  to  account  for  the 

u:e  of  diffwr.t  Hganthxic  bases,  are  -.iproxlattely  5,  S/2,  ana  2  re  spec- 

\  t/2  2 

lively.  Her-*  as  expected,  the  tliaes  in:r«aee  *s  W ,  k  and  k  respective¬ 
ly.  for  the  ncn-llr.esr  Ci'es  r.lout  10  tiojs  as  such  uschine  tiae  «*i  required 
ee  for  tl*  linear  cases. 

The  rusher  tf  iters*. tens  required  using  the  Crank-Wlcolscei  nathod  seer* l 
to  be  nearly  Independent  of  whether  the  differential  esuaticn  was  linear  cr 
oosv-llnear  or  wither  BT  or  SD  boundary  condition*  »«re  used.  The  following 
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.  « 
l 


UH,  U»  ..W.  cf  «  u..d  for  ..ch  «  »d  «...  Wi-U 


iter.tiona  required  yr  r». 


Io.  of  Iteration 


1*025 

1*069 

1.144 

1.375 


S.6.  liaUilSL  IfiJKlMioaa.  Con.id.ring  the  «naly.i.  «  the  r..uli.  a. 
a  al.clc  ve  an  led  to  th.  fcllcnring  tentati*.  conolu.icn.* 

1.  In  the  linear  case,  tr*  error  caused  by  th.  «.  cf  ai*  of  the  differ¬ 
ed  equation,  rat.ar  than  the  differential  e.uati.nc  eari.a  a.  h*.  to  th. 

non-linear  case  «As  ia  true  for  the  and  CH  r«,thod,  but  for  the  CIW(l)  and 
■stboda 

CJCJ(Il)Ath*«  error  rsries  as  h. 

*.  Ti>e  «  cf  Option  I  and  th.  cordiV.cn  (6.1.H)  us*  at  x/k  -  1  intro¬ 
duce  error,  of  order  h.  Cn  the  other  hand,  the  ua.  of  Opti*  II,  M  it  can 
te  ueed,  and  OV.lcru  III  and  III'  introduces  error,  of  order  h2. 

further  *crk  ir .odd  include  tl*  following! 

iht*  trv  CWJ  with  »  fuced  r.ucber 

1.  In  the  c.ir-lin.ur  case*,  n»«  and  w»*  ‘*5  WM* 

of  it. ret  ions  independent  of  h. 

,  HO  oat  .an  o,.u«  in’  f<»  e»  1,rl  b“"u,T 
,2.1.1S)  -  12.!. 17)  t.  r.pr,»nt  th.  rlsht  Und  h~*l*W 

s.7.  alts4-  *  *•"  “*  “*  f,"iM  *Ush  "" 

mim  a.  a»  -t»rf.  Th.  m.it.  •»>■*  «•"  <*■**•*  *tat 

tu  wtnai  .~w  <*-.  u  ..  ..  a.  «  «  cm  -u-.  *“* 

tn«y  b.  -'.ru-r  it.  »W»1*  <*  —  “*•*«  **• 
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List  cf  itbreTiatlota 

1.  f Pi3ffiU^IJfcL5ft^L£Pl  Bcuxierr  Condition* 

The  function  f(x,t,v) 


JbfffcSl 

^  tl-D 

£SH^XZ  ggRgiUSPg 

BVL 

0 

Sea  (3.1.1)  eith  g(x)*l 

fvh 

au 

Sm  (3.2.1). 

NHL 

tea  (2.1.2). 

NUN 

,-Va 

See  (2.1.2). 

for  both  KDL  and  NON  tlw  initial  value?  used  are  ^ivan  in  Table  1. 

2.  QlffrperKe  fetation? 

ID  loroard  Difference 

>B(|)  uees  r  «  | 
rD(i)  uses  r  -  J 

C*  Crank-Nlcoltcn  (iteratiro) 

CN(I)  uees  r  •  J.  1,  2,  4  eith  M  -  8,  16,  32,  64.  reepectlv*.^ . 

GSB  Crark-Niccleon  (noo-lisral in) 

CW(I)  uses  r  -  i,  1,  2,  4  *iUi  11-8,  18,  Z2,  64  respectleely. 
001(11)  uses  r  -  1,  Z,  i,  8  Mbh  a  -  8,  16,  12,  64  re?pectl»»ljr. 

Ok>  Du>  ort-f  rar.ke  1 

3-  %£D!iSllV'li«rA-9L  -  (» C4-  NO  problem  only) 


Obtlcn 

fcild-liPQ 

X 

(2.1.5) 

II 

(2.1.7) 

III 

(2.1.11] 

III' 

(2.1.13] 

TABU  I 
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MUfflJSS  Or  aurtxwc*  OtUTia  » 


mm 

asiafr 

4  aSia 


3lfil 


i3i|ll 


§s§Ia! 


3S||||  "SsMl 


WHOOp Q 

P  PVPtg, 

Si 

X3|?fc^ 


58*^.o~’ 

1"I39| 

3e8?3?’ 

lO  H  f  9  r?H  r4 

H  I  H  •  H  H 

«0  «  ^  «OMO^ 

•  iC-  HW^ 
«  45  *r  i  >n  n 
•~it*  •*  VH 

a ■ « sda 


ssgSlt 


Si 

-t  *4^t 

1  'ar* 


J 

▼  j 

« 

T-iass?  T 

¥  tfll 
•• 

f _ 

issis?  , 

M  NM<^| 

Mil 

S?H8  - 


!  WI 

«  8  P  Pid  fi  8 

i  H*& 

i  *» 


'  5  slgafcjfe 

!  H  m 

»,-S  *:J— i— r- 


f| 

H233J 

aa9>t 

I  r** 


|^x^r 


--r-T-t 


- • - 1 - _ — — 

aaRsrf?' 

O  «  m.  H?trt 
fC  M\h  f» 

™ I — ^-*z; 

383637' 

i 

WU| 

J4S1 

BOI*  S Kffll 

gaa^as 

t  i  ■-Jf'  r> 

’I  '*..•*•** 


~m. 


ISciS? 


pslli 

*?3wa 

-  35  ^*»Cj 


rajaa  * 

TgsasE' 

r&M,  - 


laifH 


i”  i 


rm  1 1 

i  7 

ai&3£  v  ! 


*  •  #i  *  i  V 11  • 


»***•* | 
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Appendix 

Fomulaa  Vwd  for  Ortrec.  Proem 


Notation* 


a  -  ms  dlacusaed  In  Section  S 
if  J 

(n) 

u.  .  -  value  of  function  for  n-th  iteration 
'tJ 


Bracketed  tern  included  only  for  the  non-linear  caee 
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KB  -  W 
« 


,  ♦  U  -  2u  )  ♦  tt  ♦  I  lrf  i*  M-i;  3>< 

Uyl  •  HVl,J  4  Vl»i  2°^  1.1  ^  j* 


'•**  “  jae  o 

• 

ND 

n  * 

i,J*l 

r  -  Vni,j1 

r(«  ♦  u  “  2ti  .)  ♦  o.  4  4  1  to  I  * 

'  i-l,J  i*l,J  M  L.  -* 

i*i«it-r,  >o 

Op  I 

.  »*"«  * 

*0,J41  1  ♦  Hh 

J2*0 

op  n 

«  •  u  ♦  Kh(a  -  J  ♦Fka 

ofJ*i  u  *  ^  -» 

J*0 

Op  III* 

r  •  J  only 

•i  ,.i ♦“»„♦»{> j'C* 

!]  ... 

°»Jn  1  ♦  Kh  ♦ 

All  Optional 
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cm  -  BT 


d^.«a  f“ 

\  “  4  °2,j  -  2ul,j>  4  ir  * 1,1  *  ?*?  4[li?kT*  'Ll  J*° 

dl  •  zfeW-lJ  4  Vl,J  -  2ui,j}  4  A  *i.J  4|lt  J*0 


dU_l  '  2*2r(V2,j  4  Vj  ■  2Vl,J5  4  1*T  Vl,J  4  Z*r 
^  *  0 


b  "0 

0 


b1  “  2*2r 

-  r/2*2r 

b‘ ' » -fefejv. 

b  •  0 

U-l 


2*  i*U -2 


qi  "dl 


dl  *  (z+2r, 


(z<ir)\-l 


24  1*1 M 


Vl,  J4l  "  Vl 

'S’  biVi,j«i' 


IS  IS 


dc 

\  ‘  *  \3  ‘  2Ul,iJ  +  ^  2^  Vi*  tE* 

di  -  zfe^j  4  3*  j  ■  *4$  *  &  Vi  ^  k*  ,3J 

S.'0  . 


b  -  0 
0 


-  r/2+2r 


-  0 


*»  jj 


Vj*l  ’  ^ 

tti.  i*i  *  *i  “  Vi*, i* 
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OFT 

ft  A  « 

Compute  2  rows  by  iD  technique  using  r',  where  r  *  r'*2  |  O^r'^). 

Let  2®  th  row  by  KD  be  J  •  1  (first  row)  for  D»K.  Them 

tti,J*l  “  I^Vl, j  4  Vl,J  "  2ui,J-l*  4  *1,J-1  4  jj?***^/ 

UKK-l)  jjfr 


4 

where  e 


BT 


MD 


The  boundary  conditions  are  the  same  as  for  the  ID  method,  both  BT  and  MO. 
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